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A COMPUTER PROGRAM FOR SYSTEMATICALLY ANALYZING 


FREE -FLIGHT DATA TO DETERMINE THE AERODYNAMICS 
OF AXISYMMETRIC BODIES 
By Gerald N. Malcolm and Gary T. Chapman 
Ames Research Center 


SUMMARY 


The computer program analyzes free-flight motions of axisymmetric bodies 
to deduce coefficients of drag, lift, and static and dynamic stability from 
any set of free-flight data complete enough to define the trajectory and angle 
history. Nonlinear behavior can be accurately assessed from the results of 
the data-reduction program, and available methods are discussed. To demon- 
strate the effectiveness of the data-reduction method, free-flight tests of 
the AGARD standard hypersonic ballistic correlation model HB-2 were conducted. 
The results for lift and static and dynamic stability at a Mach number of 2 
agreed well with conventional wind-tunnel results. Significant differences 
were found in the drag data because of different base pressures believed due 
to sting effects in the wind-tunnel tests. 


INTRODUCTION 


Free-flight techniques have been used at Ames Research Center, in both 
ballistic ranges and counterflow facilities, to obtain aerodynamic characteris 
tics of many types of configurations. In past years, raw data were reduced 
principally for drag and stability and only occasionally for lift. The 
procedures involved several computer programs and a substantial amount of 
manual calculation. 

Recently, a Fortran IV computer program was developed to systematize the 
data-reduction procedure and eliminate all manual calculations. The program 
produces in one operation a set of quasi -linear values 1 for the aerodynamic 
coefficients of drag, lift, and static and dynamic stability. Improved proce- 
dures have been incorporated to calculate drag coefficients, and the method of 
determining lift and stability coefficients has been modified to reduce errors 

To demonstrate the effectiveness of the data-reduction method, it is 
applied to free-flight data from tests of the AGARD standard hypersonic 


1 The values are termed "quasi -linear" since the equations of motion that 
are solved explicitly assume the static forces and moments are either constant 
or vary linearly with angle of attack and dynamic forces and moments vary 
linearly with angular rate. 



ballistic correlation model HB-2, and the results are compared with 
conventional wind-tunnel results (refs. 1-3)* 
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moment of inertia about the axis of symmetry 

moment of inertia about a transverse axis through the center of 
gravity 

constants in equation (4) 
mass of model 
Mach number 

roll rate about axis of symmetry of model (positive clockwise 
looking from rear) 

pitching velocity 

free -stream dynamic pressure 

Reynolds number based on d and free -stream conditions 
true distance along the flight path 
flight time 

velocity along flight path 

range coordinate in direction of flight path 

horizontal range coordinate perpendicular to x and z axes 

vertical range coordinate perpendicular to x and y axes 

angle of attack (in the vertical plane) 

average value of maximum-angle envelope 

average value of minimum-angle envelope 

/i r x 

root -mean-square angle of attack, / — J ocr dx 


“r 

a R 

P 

T ]l,2 
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small -angle approximation to resultant angle of attack, \l a 2 + P 2 

resultant angle of attack, tan -1 \J tan 2 a + tan 2 p 

angle of sideslip (in the horizontal plane) 

damping exponents in equation (4) 

angular displacement measured in the xz plane 
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i 


dynamic -stability parameter, Cp - Ol + 



p free -stream air density 

a radius of gyration about a transverse axis through the center of 

gravity of a model 

$ angle of rotation of range coordinate axes about the x axis 

\|r angular displacement measured in the xy plane 

toi 2 rates of rotation of vectors that describe the model oscillatory 

motion in equation (4) 


DATA -REDUCTION SYSTEM 


To obtain aerodynamic coefficients from a free-flight test, the basic 
data necessary for analysis are position of the center of gravity, angular 
orientation, and time at various points along the trajectory. The earth-fixed 
orthogonal coordinate system in ■which these data are desired is shown in fig- 
ure 1. Other requirements are a knowledge of the mass, moments of inertia, 
center-of -gravity location, and geometric characteristics of the model, 

together with the physical properties 



x - axis in direction of flight 
y- horizontal axis (positive right) 
z- vertical axis (positive down) 
x B - axis through of body 
9- projection of angle x 8 o x in x-z plane 
f - projection of angle x B ox in x-y plane 
p- roll rate about x B axis 

Figure 1.- Earth-fixed orthogonal coordinate 
system . 


of the test 
is fired. 


gas into which the model 


The data-reduction program 
includes the following routines: 

1. Convert raw measurements into 
earth-fixed orthogonal coordinates 

2. Determine drag coefficient 
and gravity corrections 

3. Determine static and dynamic- 
stability coefficients 


4. 

system 


Rotate orthogonal coordinate 


5. Determine lift coefficients 

6. Correct measured angles for 
flight -path swerve 


A description of each of these routines and their dependence on one another 
will be given. 
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Obtaining Measurements in Earth-Fixed 
Orthogonal Coordinate System 

The raw measurements obtained from a particular facility must be proc- 
essed to whatever extent necessary to express them in a set of earth-fixed 
orthogonal coordinates . This procedure depends on the fiduciary system of the 
test facility from which the data were obtained and would be bypassed if the 
data were in the proper coordinates initially. 


Drag and Corrections for Gravity 

The drag coefficient is obtained directly from the flight time and dis- 
tance measurements by using a high-precision least -squares curve -fitting pro- 
cedure (see appendix A). The equation relating time and distance can be 
written (ref. 4): 


1 e KC ^ 

t " to “ V d KC d + V 0 KCj) (1) 

where t = t 0 and V = V Q at x = 0 and K = pA/2m. The parameters C D> V Q , 
and t Q are determined to give the best fit to the experimental values of x 
and t. Once these coefficients have been determined, the position and angular 
measurements in the vertical plane are corrected to account for the influence 
of gravity. A procedure is also available for eliminating erroneous time 
measurements (discussed in appendix A). 

To this point, the drag coefficient has been assumed to be constant for a 
given flight (or segment of flight) and not a function of angle of attack. 

What has been determined is really an effective drag coefficient over the 
angle -of -attack range encountered in a given flight. From several flights at 
different amplitudes, the dependence of drag on angle of attack can be deter- 
mined. For most bodies of revolution, the drag as a function of angle of 
attack can be expressed as 


°D * C D 0 + < 2 > 

where 0 d o is the zero -angle drag coefficient and ap is the resultant angle 
of attack. It can be shown (ref. 5) that the effective drag coefficient we 
have determined is 


^D 0 ff ~ ^D Q + (3) 

where 
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From equations (2) and ( 3 ) , then, Cp^^ is the value of Cd that occurs at a 
resultant angle of attack equal to the root -mean-square resultant angle of 
attack of the flight in question. Therefore, for each flight or flight seg- 
ment, one obtains a value for Cn e ff and arms* If these points are plotted as 
CDeff versus a§m s and a straight line is fitted to them by least squares, 
values for Cp 0 and Ci are obtained , and one then has an expression for Cp 
versus angle of attack. 

If the drag coefficient cannot be expressed as a quadratic function of 
angle of attack but is, instead, of the form Cp = Cp 0 + C 2 |aR| n , then the 
angle of attack at which the effective drag coefficient should be plotted is 



° D eff = Cd o + ° 2 ^ (ref * 6) 

The best value of n is found by a trial-and -error process. 


Static and Dynamic Stability 

Stability coefficients are determined from the oscillatory history of the 
model. The well-known tricyclic equation derived by Nicolaides (ref. 7)j mod- 
ified to use distance rather than time as the independent variable, is used to 
analyze the angular motion (see appendix B) . Distance is used as the indepen- 
dent variable because it essentially eliminates (to terms of second order) the 
dependence of the calculations on velocity variation and hence is more conve- 
nient. The important assumptions of this method are linear aerodynamics; 
small angles; constant roll rate; small asymmetries, mass and configurational; 
and small velocity change. The modified equation is 

3 + ia = Kl e (,ll+i “ l)x + K 2 e ( '' 2 - 1 “ 2)x + Kse 1 "* (M 

where 

Ki = bi + iai 

Kp = b 2 + ia 2 
K3 — b3 + ia3 
P = (“1 - W2)(ly/lx) 

A least-squares procedure using differential corrections (ref. 8) is used to 
fit equation (4) to the experimental data, a and |3. The ten constants bi, aj., 
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b 2 , &2) b 3 , q.2,) tu* <*>i, t) 2 , and co 2 are determined from the fit. As shown in 
appendix B, when the Magnus moments can he ignored (generally, the case for 
most ballistic-range tests unless the model is deliberately spun) , the con- 
stants r)!, r\ 2 , a)!, and w 2 are related as follows: 


T)! 4- (pA/2m)Cr 
0)2 t ]2 + ( pA/ 2m) C 

where CL a is the lift-curve slope. Because o) x and o) 2 are easier to deter- 
mine accurately, rj 2 can be written as a function of w 1; and oo 2 , thereby 

reducing the number of unknown coefficients to nine (0L a must be known and 
its determination will be discussed later). If p (the roll rate of the model 
about its own axis, assumed constant) can be measured in flight, u>i and w 2 
are not independent of one another and the number of coefficients is reduced 
to eight. The aerodynamic parameters of static and dynamic stability are 
related to the determined coefficients as follows: The quasi-linear pitching- 
moment -curve slope, , is related to and as 


-2I V 

■ ysa < 6) 

and the dynamic-stability parameter, £, where £ = Cp - Cx, a + (d/cr) 2 (C + C^) 
(ref. 9) is related to tj x and r \ 2 as 


Til + T)g 
pA/2m 


(7) 


The remaining six coefficients, a^ and bj_, along with r\ x and r} 2 , describe the 
envelope of oscillatory motion for the particular flight being analyzed. 

The assumption of linear aerodynamics, that is, static forces and moments 
that vary linearly with angle of attack and dynamic forces and moments that 
vary linearly with angular rates , does not prevent the use of the method for 
bodies with nonlinear stability coefficients. In such cases, the method is 
used to reduce data from several flights or portions of flights at different 
amplitudes. These quasi-linear values for various angle -of -attack amplitudes 
are then used to obtain the desired coefficients as functions of angle of 
attack. The method of obtaining nonlinear pitching moments from quasi-linear 
data is derived in reference 10 and illustrated in some detail in reference 11. 
Basically, for a pitching -moment equation of the form 


-Cm = p Q a + pi a 2 + p 2 a 3 


the equation for the quasi-linear 
, can be written 


+ p 3 a 4 + p 4 a 5 + psa 6 + Psu 7 * * * (8) 

value of the pitching -moment -curve slope, 
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( 9 ) 


where 


“m + “min /°m + “min'' 2 


h = 


^“m + “minj 


(cm “ ctmin ) 2 


c = °m + “min 


“m + “min f ° m + c 4in'' 7/2 


d = 


(“m " °min) 2 


e = 7< + 10c§4 in + lo^ n 


“m +a mm /“m + “min' 9/2 


f = 


(4 - “min) 2 


g 3°^ + 5a m a m j_ n + 5a m a m ^ n + 3 a -mm 


and c% and Om -fn are the average maximum and minimum resultant angles, 
respectively., in the a-0 plane (sketch (a)). 



£ 


Sketch (a) 


Values for C™ , and from several 

flights are fitted by a least -squares procedure 
using as many terms of equation ( 9 ) as desired (a 
computer program has been written to accomplish 
this systematically (ref. ll)), and the resulting 
coefficients p 0 j Pij P 2 j • • •; etc., are deter- 
mined. These coefficients then produce an expres- 
sion for C m versus a. A method is also available 
to determine nonlinear results for the dynamic - 
stability parameter (ref. 12). It has been the 
author's experience, however, that the quasi - 
linear values for dynamic stability cannot be 
determined accurately enough to define the 
nonlinear contribution. 
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Rotation of Orthogonal Coordinate System 

One of the assumptions necessary in the derivation of equation (4) is 
that angular displacements are small and that the resultant angle of attack is 
simply the square root of the sum of the squares of the pitch and yaw angles: 


ct£> = \l a 2 + (3 2 


( 10 ) 


The exact expression for the resultant angle of attack is 


= tan -1 v tan 2 a + tan 2 (3 


(id 


These two equations agree closely for small angles. They agree exactly if a 
or 3 equals zero. Thus., the error introduced hy equation (l0) can he reduced 
for nearly planar motions hy rotating the coordinate axes through an angle 0 
so that the motion occurs near either -the pitch plane or the yaw plane, thus 
keeping the angles in the other plane small (see fig. 2). The procedure for 
rotating the data is described in detail in appendix D. 


8 — 


6 — 


4 — 


2 — 


a, deg 

0 


-2 — 


-4 — 


-6 — 


- 8 *— 



deg 


The error due to use of equa- 
tion (l0) for large -amplitude motions 
is most pronounced in the dynamic - 
stability parameter | . Figure 3 
(reproduced from ref. 13 ) demonstrates 
the magnitude of errors encountered for 
a typical case. We consider a test in 
a range having 11 data stations at 



Figure 2.- Rotation of coordinate axes. 


Figure 3*- Effect of resultant angle assumption 
on dynamic stability parameter. 
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4-foot intervals with test conditions such that pA/2m = 0.0006/ft and the 
wavelength of oscillatory motion is 26 feet. We consider various values of 
arms? and £. Figure 3 shows the induced error (| - f; exact) for $ = 45° 
(the worst case) as a function of £ for three values of cy ms . Note the 
near linear dependence of the induced error on |, and also the strong 
influence of ar ms • 


Lift 

The lift-curve slope, CL a , is determined from swerve measurements (with 
effects of gravity removed) in conjunction with the oscillatory motion of the 
model. A modified form of Nicolaides' equation is fitted by the method of 
least squares to the experimental displacement data z and y. The equation 
used is (see appendix C for derivation) 


y+iz = - |^C L 


PX PX 

'o 


( (3 + ia)dx dx + (d)(CL q + CL a 


>x 


+ ( c yo + iC L 0 ) 


1 + ipx - e^l x ' 


ipx^| 


u 

+ (yA + iz')x + (y Q + iz G ) 


( 0 + ia)dx - ( Po + iao)x 

( 12 ) 


The integrals are obtained in closed form by integrating the tricyclic equa- 
tion for 0 + icc using the coefficients solved for in the stability routine. 
Therefore, the constants CL a > (CLq^ + CL<x) , Cy Q , CL 0 , yo, zA, and z 0 
appear in a linear fashion, and a straightforward least -squares technique can 
be used to determine them. 

Although so far most data have been reduced only for linear lift coeffi- 
cients, it is fairly straightforward to include suspected nonlinearities. For 
example, if the lift curve can be approximated by 


CL = C Lc a + Cl 2 cx 3 


one additional term should be added to the right side of equation (12), namely 
(ref. 13), 


- c l 2 f o f o Ip + lct ! 3 ^ dx 


The additional constant also appears linearly and can he solved for as before. 


Corrections to Measured Angles for Flight -Path Curvature 

Because the original angular measurements are made with respect to 
earth-fixed axes and not to the actual flight path, a correction to the angular 
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data may be required for the swerving of the trajectory as derived from the 
first pass through the lift routine. Figure k shows the relationship of the 



y 



(b) Horizontal plane. 

Figure 4.- Relationship of swerve corrections and 
original measured angles in vertical and 
horizontal planes . 


swerve trajectory, flight -path angle, and measured angles in orthogonal planes. 
A detailed description of angular corrections is outlined in appendix E. Once 
this correction is made, the data are again cycled through the stability rou- 
tine with the corrected measured angles, and this process of recycling 
continues until the swerve corrections to the angles remain essentially 
constant . 


Program Summary 

When the individual routines are completed, a final tabulation of the 
calculated angles and displacements along with the computation of the root- 
mean -square angle of attack, a^g (necessary for representing the drag 
coefficients), is made for a predetermined array of x values that can be 
used for plotting and reference. 


Ss.- 


ll 



A typical procedure for reducing a free -flight run can he summarized as 
follows (see fig. 5 for flow chart): 


Raw data 

1. - Film measurements 

2. - Model characteristics 
3- Test conditions 



1- Obtain measurements in earth - 
fixed; orthogonal coordinates. 

2. Fit the x-t data with equa- 
tion (l) and determine the drag coeffi- 
cient. Check for bad time data with a 
smoothing routine and recalculate the 
drag coefficient if necessary. With 
time now given as a function of x, 
compute the vertical drop of the model 
due to gravity and modify the z and 0 
measurements to eliminate its effect. 

3 . Use the Q, 9 and x data in 
the stability program to determine the 
orientation angle ($) of the motion 
with respect to the xy plane. 

4. Rotate the Q, z> and y 
values through the angle $ so that 
the motion occurs as much as possible 
in the xy plane. 


Figure 5 .- Flow chart of data reduction procedure. 5* Rerun the 6, l]/"; and X values 

through the stability routine to obtain 
the first set of calculations for the constants in equation (4). 

6 . Use the y, z, and x data in the lift program along with the computed 
coefficients from the stability routine; and determine the constants in 
equation ( 12 ). 


7* With this set of constants from the lift program; one now has expres- 
sions for z versus x and y versus X; which can be differentiated with respect 
to x to obtain corrections to 0 and *i|r and thus to estimate a and ( 3 . 

8 . With this new set of angles; the stability program is called again 
and steps 5 through 7 are repeated until the computed angle corrections are 
essentially constant. One then has the best fit to the actual oscillatory 
history with respect to the flight path. 

9* At this point the a; |3; z } and y data may be checked for data points 
that appear to be incorrect (appendix F) . If any are found; the values are 
replaced by calculated values and the process of iteration between stability 
and lift routines is repeated as before. 

10. All measured and calculated values for a; {3; Z; and y are finally 
rerotated to the original axis system of reference. 
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The program offers a degree of flexibility to the user; that is, some of 
the procedures can be done independently of the others, such as the drag rou- 
tine. The rotation and swerve correction routines are optional as well as the 
smoothing routines for eliminating bad data. The lift routine can be run inde- 
pendently if one provides the necessary coefficients from the stability analy- 
sis of the angular motion. The computing time to reduce one run or run se gm ent 
is typically 0.5 min on an IBM 7094. Other minor but very convenient (and 
sometimes necessary) features can be utilized with the program. Additional 
details of the program are available on request. (Write to Ames Research 
Center, NASA, Code SVH, Moffett Field, California 9^-035) ■ 




Solid model 

Note; All dimensions in inches 



(b) Two types of models. 


COMPARISON OF FREE -FLIGHT 
AND WIND-TUNNEL 
RESULTS 


Tests were conducted at 
a Mach number of 2.0 and a 
Reynolds number of 1.7 million 
in the Ames Pressurized Bal- 
listic Range with the AGARD 
standard hypersonic ballistic 
correlation model HB-2 (see 
fig. 6). Where possible, the 
results are compared with 
results obtained in a conven- 
tional wind tunnel (refs. 1 
and 2) and with results from 
ballistic -range tests 
(ref. 14) . 

The two types of models 
used in the Ames tests are 
shown in figure 6(b). The 
first set of models was solid 
in construction and was made 
of phosphor bronze and a 
polycarbonate -type plastic 
(Lexan). A second set of 
models of the same materials 
was constructed to have a 
smaller mass and transverse 
moment of inertia in an unsuc- 
cessful attempt at obtaining 
high angles of attack. This 
set of models had a hollowed 
flare and cylinder section as 
shown. Two pins were mounted 


Figure 6.- Sketch of model. 
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in the base of both sets of models to permit the measurement of roll in flight; 
therefore , the roll rate was a fixed input to the data -reduction program. 

The data used for comparison with wind-tunnel results will he that 
obtained from the solid models. (The data for all the tests are given in 
table I.) To reduce the effects of aerodynamic nonlinearities with angle of 
attack, the data were analyzed with all flights split into overlapping segments 
of approximately l-l/2 cycles of motion with three amplitude peaks. Figure 7 
shows a typical set of position, time, and angle data for the present tests. 
Figure 7(a) shows the flight time as a function of distance and the resulting 
velocity decrease over the length of the range. Figures 7(b) and (c) show the 
displacement and angular measurements ( z and 9 were corrected for gravity) . 




Figure 7.- Typical flight data. 
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Figure 7.- Concluded. 


Estimate of Error 


A method for computing the effects of random errors on the various aero- 
dynamic parameters measured in a ballistic range is outlined in reference 13 * 
With the equations contained therein, estimates of the standard deviation in 
the respective parameters due to experimental error were made assuming a uni- 
form error distribution with maximum random errors in the measured data of 


y,x) = ±0.001 ft. 

and (t) 

= ±1X10 6 sec. 


a = 2° 

a = 6° 

sd(c d ) 

± 0 . 005 

±0.005 

SD(CiJ 

±•3 

±.l 

SDtc^j) 

±.009 

±.003 

SD(g) 

+ .6 

±.2 
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Drag 



a, deg 


Figure 8-- Variation of drag coefficient with 
angle of attack- 

ill the wind-tunnel "base pressure meas 


The total drag coefficient as a 
function of angle of attack is shown in 
figure 8* 2 3 Each data point was obtained 
by analyzing a flight (or portion of 
flight) as already outlined. These 
data points were then fitted by a qua- 
dratic equation in angle of attack by 
the method of least squares. One imme- 
diately notices the significant differ- 
ence between the present ballistic -range 
results and AEDC wind-tunnel results. 

The total drag for the wind-tunnel 
results was obtained by combining the 
contributions of forebody and base drag 
measurements given in reference 1. It 
was felt that the discrepancy could be 
*ements since sting effects on these mea- 



2 Most of the scatter in the data points can be attributed to slight dif- 
ferences in Mach number. Xh table I, the drag coefficient increases slightly 
with a decrease in Mach number for a given range of cym^ . 

3 The intense shock pattern shown on the last half inch of the model is 
the area where the sabot gripped the model. 
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of the wake, the Mach number of the flow along the wake can be determined and, 
consequently, the ratio of static-to-total pressure in the wake (and on the 
base) is known. The total pressure behind the normal shock at the nose is 
assumed to be constant everywhere on the body and, as a result, the ratio of 
static pressure at the base to free-stream pressure is known. The base drag 
can then be calculated. The calculated values from all four flights were aver- 
aged and added to the AEDC forebody drag; the shaded curve shown in figure 8 
is the result. The total deviation from the average base drag is indicated by 
the shaded error band. The agreement with the free -flight total drag is now 
excellent . 


Lift 


The lift -curve slope, as deduced from the swerving motion of the model, 
is shown in figure 10. These data indicate that the lift coefficient is lin- 


4 

/— AEDC 
-L- O 

- - r— C 0 

C L<t . per 

radian 5 

— fc) 
^^AEDC with Ames 
— base drag 

w W O 

U CLa = 3. 40 

2 

1 l 

1 1 


Otm, deg 

Figure 10.- Lift- curve slope. 


ear with angle of attack at least to 7°* 
The wind-tunnel result for lift -curve 
slope, as deduced from the normal and 
axial force measurements, is shown, 
together with a value obtained using 
the Ames deduced base pressure in the 
axial force contribution. The agree- 
ment between free -flight and wind-tunnel 
results is very good. 



Normal Force 

A curve of normal force versus 
angle of attack was calculated from the 
measured lift and drag data where 
Cl = 3-400. and Cd = 1-239 + 0.00249a 2 
and ON - CL cos a + CD sin a. The 
results are plotted in figure 11 and 
are compared with AEDC wind-tunnel data 
and with CAEDE ballistic-range data 
(ref. Ik ) . 


Figure 11.- Variation of normal- force coefficient 
with angle of attack. 


Pitching Moment 


The nonlinear pitching -moment - 

coefficient curve was calculated with the method of reference 9> with linear 
and cubic terms in angle of attack chosen as most representative of the data. 
(Higher order polynomials were examined, but were nearly identical to the 
simpler linear-cubic representation.) One flight at low amplitude (am = 1.8°) 
revealed a significant decrease in static stability (see table I, run 1156), 
which is reflected in the lower initial slope of the pitching -moment curve in 
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Figure 12.- Variation of pitching -moment 
coefficient with angle of attack. 
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(b) Damping in pitch. 



Figure 13- - Variation of dynamic aerodynamic 
coefficients with angle of attack. 


figure 12 . Curves from AEDC and CARDE 
are also shown and indicate slightly 
higher values for the pitching moment. 
Small differences in boundary-layer 
conditions and base pressure could 
contribute to the discrepancy. 


Dynamic Stability 

The dynamic -stability parameter 
for nonpowered flight at constant alti- 
tude, | = Cj) - C]^ + (d/0) 2 (C m ^ + Cjq^) , 

is shown in figure 13(a) and the result- 
ing values for (Cm^ + C^) in fig- 
ure 13(b). The value of (C^ + C^) 

obtained in conventional wind-tunnel 
tests (ref. 2) is slightly lower than 
free -flight results. The dynamic con- 
tribution to lift due to pitching and 
plunging (Cl^ + Cp^) , which can also be 
considered as the contribution to 
damping in pitch due to lift (other 
major contributions would be center-of- 
pressure movement or changes in drag 
due to angular rates), is shown in fig- 
ure 13(c). To the author's knowledge, 
this is the first fairly consistent set 
of data obtained for ( Cl^ + QL<x) in a 
ballistic range. 


Effects of Hollow-Base Models 

While, as noted earlier, the 
attempt to obtain higher angles of 
attack by reducing the mass and trans- 
verse moment of inertia of the models 
was not successful, the data from this 
set of tests revealed some interesting 
results . Figure 14 shows the drag 
coefficient for both solid and hollow- 
base models. Although some of the 
scatter in the data for each model type 
can be attributed to slight Mach number 
effects, the hollow -base models had a 
slightly lower overall drag coefficient 
than the solid models. The reason for 
the lower drag is not known but may be 
due to subtle boundary-layer differences 
or effects of the hollow base or both. 


18 



1.6 


-D 1.2 


.8 




•5 




• Solid models 
□ Hollow models 


Data for the other aerodynamic param- 
eters for the hollow models were very 
similar to data for the solid models 
except that there seemed to he more 
scatter. It is possible that the open- 
base creates a "sloshing" effect of the 
air movement in that vicinity and may 
introduce slight differences in the 
forces and moments. It should be 
pointed out that an open-base model was 
used in the CAKDE tests , but the enve- 
lope of oscillatory motion for the 
tests never reached below approximately 
7° and therefore might not be sensitive to these effects. Model base geometry 
should be carefully considered in the design of models for tests at subsonic 
and moderate supersonic Mach numbers at which base pressure is comparable to 
the free-stream pressure and the base forces are comparable to forebody forces. 


4 

a, deg 


Figure 1^.- Drag measurements for solid and 
hollow-base models . 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif. 9^035, May 16, 1968 
124-07-02-35-00-21 
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APPENDIX A 


DRAG COEFFICIENT AND GRAVITY CORRECTION SUBROUTINE 


The purpose of this subroutine is twofold. First, it obtains a least- 
squares fit of the drag equation to a set of measurements of distance flown as 
a function of time to deduce a drag coefficient. Second, the drag equation is 
then used to calculate corrections to the displacement and angle measurements 
in the vertical plane to obtain a zero gravity trajectory. 


Drag 


The equation of motion along the flight path is 


m 



_ 


2 °D\dt/ 


dsY 


(Al) 


where D is the drag and s is the distance traveled by the model center of 
gravity. The present analysis is based on the following assumptions: 

1. The distance s can be replaced by x, the distance along the earth- 
fixed axis in the direction of flight. This assumption is very good since the 
lateral motion is small compared to the motion in the x direction. 

2. The air density, p, the model mass, m, the reference area, A, and the 
drag coefficient, Cp, are constant. 1 

With these assumptions and the initial conditions t = t 0 and V = V Q at 
x = 0 the solution to equation (Al) is 


1 e KC P* 

° " VoKC d + VoKC D 


(A2) 


where K = pA/2m and the subscript o denotes initial conditions. 

To obtain a least -squares fit to the x and t data, we write 
equation (A2) as 


t = Ae CX + B (A3) 

1 In most tests the velocity loss during flight is small; hence, changes 
in drag coefficient because of diminishing Reynolds number or Mach number may 
usually be neglected. (This may not be true in transonic tests or, to a small 
degree, in the present tests.) In most cases, however, the drag coefficient 
is a strong function of angle of attack. This point is considered in the text 
and a method for detentlining the variation with angle of attack is given. 
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where 


VoKC D 


B 


t o 


1 

v q kc d 


c = kc d 

We will now set up the least -squares procedure. The sum of the squares 
of the residuals,, SSR, is written 


SSR 


N 

^ (^expq 
i=i 


tq) 2 W ± 


(A4) 


where t e xpq is the experimentally determined time at the ith position, tq 
is the calculated time, Wq is a weighting factor, and the summation extends 
over all positions N. (The normal mode is to operate with a Wq of 1; how- 
ever, one can use 0 as a mechanism for rejecting erroneous data points.) 
Substituting equation (A3) into (A^-) yields 


SSR 


N 

^ (texpi " Ae Cxi - B) 2 W ± 
i=i 


(A5) 


The normal least -squares procedure would be to form the partial derivatives of 
equation (A5) with respect to A, B, and C and set them equal to zero, and 
solve for A, B, and C. This works as long as the partial derivatives are 
linear in the unknown coefficients, but it does not work here because C 
appears in an exponent. Hence, we will proceed as follows: Form the partials 

with respect to A and B and set them equal to zero , yielding 

N 

^ WpCtexpi “ Ae Cxi - B)e Cxi = 0 (A6) 

i=i 


N 

^ Wi(texpi " AeCXi - B) = 0 (AT) 

i=i 

Using equations (A6) and (A7) eliminate A and B from equation (A5) to yield 


dSSR = 
Sa 


5ssr 
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SSR = 



We now have an equation in C only. At this point, a differential correction 
procedure (ref. 13) could be used, but the partial derivative with respect to 
C would be quite complicated. Therefore, a numerical procedure is used to 
minimize the SSR as follows: An approximate solution for Cp and hence C 

is found (to be discussed in the next paragraph), and equation (A8) is evalu- 
ated. Then the value of Cp is increased by e (e.g., e = 0.1 Cp) and equa- 
tion (A8) is reevaluated. If the SSR is smaller than the first value, Cp 
is increased by e again, continuing until the minimum in SSR is passed, at 
which point the size of e is reduced by one-half and Cp decreased by e. 
This process is repeated until the minimum is passed again, at which point e 
is again reduced by one -half and Cp increased by e. This process is 
repeated until Cp changes less than some prescribed amount. Normally, a 
value of 1 percent of the last value of Cp obtained is used. If, by chance, 
the first step had produced an increase in SSR, the direction would have been 
changed. Equation (A8) is evaluated to find SSR in double -precision 
arithmetic to insure accurate determination of the minimum. 

The starting solution for the above iteration procedure is found by uti- 
lizing the small -velocity-loss criterion, that is, small KCpx, and expanding 
the exponential in equation (A3) to order x 2 to obtain 

t = (A + B) + ACx + x 2 (A9) 

or 

t = A + Bx + Cx 2 (A10) 


By least -squares fitting equation (A10) to the x - t data, values of A, B, 
and C and, hence, Cp are obtained. The least-squares fitting of (A10) is 
straightforward since all three unknowns appear in a linear manner. 
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Gravity Corrections 


After convergence of the drag routine, gravity corrections can be calcu- 
lated. The new vertical position becomes 


z n± “ Zexpi 


1 ,2 
p 


(All) 


where znq is the new z at the ith position and z exp± i s th e experimen- 
tally determined z at the ith position, and the curvature of the flight 
path due to gravity is 


dz 1 _ gij 

dx I x V± 

Hence , the new angular coordinate is 


(A12) 


®exp^ 


+ 


tan" 1 


dz 

dx 


i 


(A13) 


At this point we may elect to check for erroneous data points as follows: 
The standard deviation of the time measurements is calculated from the SSR 
as 

SD(t) (Al4) 


Now the difference between the experimental and the calculated time from the 
least -squares fit is compared to SD(t). If the absolute value of this dif- 
ference exceeds a chosen multiple of SD(t) (normally taken as 2), the data 
point is deleted from the set of data. "When all data points have been checked, 
the least -squares procedure is repeated if points were deleted; if points were 
not deleted, the procedure is finished. When data points are deleted, a 
statement to this effect is printed in the output. 
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APPENDIX B 


MODIFIED TRICYCLIC EQUATIONS 


Nicolaid.es 1 differential equation of motion using trajectory-oriented 
forces rather than "body-oriented forces can he expressed in time coordinates 
as follows: 




i + v 2 



p lx pA / pd 

+ V Iy 2 m V." V UL p a 



t 


= SOOe 1 ^ dt <B1) 


where £ = 0 + ia and mb 2 = Iy. If we change to x coordinates 


Note that 


t = 


. d 


t = S£* + * + X 

5 = xg» + x 2 £" 


• ^LL dx 


x s V 


(B2) 


(B3) 


x S ^ V 2 ( -CL7 - Cd) 

where J - |/y _ i/ z | is the local flight -path angle (angle between the local 
velocity vector and the x axis of the range). (See fig. 4. ) For small 7, 
as is the case in a ballistic range where the velocity vector is nearly paral- 
lel to the x axis, the Cl 7 term can be neglected. Substituting for x 
and x into equations (B2) and (B 3 ) yields 


i = vs* 


(B2a) 


5 


v 2 


( £ 


11 



(B3a) 
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Also, 


p(t) = Vp(x) 

Substituting equations (B2a ) , (B3a), and (b 4) for (Bl) yields 


(B4) 


i 


pA . pA 

2m 6 1 2m pd 


^Pa + 


+ 


-{S 


[" 2Iy ( Cma + lpdCm Pa) + P Iy 2m 


Cm Pd " 




- P dC L Pa " iQL a 


t = Qe ipx (B5) 


vhere 


£ = C D - C La + l 1 c mn + C] 


a q ^ um a. 


Differential equation (B5) can be written as 


?" “ (9i + 92 )?' + 9x92? ~ K 3 (cp 3 - 9x) (93 - 9a)e 


93X 


which has the solution (see sketch (b)) 


ia 



5 = K 1 e' PlX + K e e' P = x + Kse’’ 3 * (B6) 


vhere 


cpi = tji + iu>x 


cp£> — ^2 “ i^2 


cp 3 = ip 


Kj. — "bi + iai 


K 2 ~ t)2 + 13-2 


K3 - 1*3 + ia3 
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Note then that 


p A 

9i + 92 = 7^ S 


i 





9 i 92 = - 


pAd 


21- 


■y 


Cm a + ipdCm Pa 


. _ tc pA 

P l y 2m 




AlSO ; 


K 3 (cp 3 - cp 1 )(cp 3 - cp 2 ) = Q 


91 + 92 ~ ’ll + ^2 + i( w i “ ^z) 


9 192 - T li T l2 + w i w 2 + i( T l2 w i ~ t 1i w 2) 


Therefore, 


T)1 + TJg = 


pA 

2m 


0)1 - w 2 = pd 


pA 

2m 


So, + W - % 


+ P 


pAd P I x pA 

1)1112 + UllJ2 = " 2lJ Sx" p d I y 2m S?, 


’a 


Xx PA 


pAd 


^ ^ - -p ^ s 0^ - pd ^ °»p a 


If one neglects Magnus forces and moments and if r iif]s « 0 )iU) 2 , then 

. _ r H + Jjg 
5 PA 

2m 


«i - w 2 = p 


*-x 


L y 


Si 


2IyO) 1 OJ 2 

pAd 


and from equation (BIO) 


c l . 


a 


£x PA 
p % 2» 


(BT) 

(B8) 

(B9) 

(BIO) 

(B7a) 

(B8a) 

(B9a) 

(BlOa) 
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However, the calculated values for t)i> t 12^ and oo 2 are not precise enough 
to obtain accurate values of CL a - A value for CL a is determined much more 
accurately from the plunging motion. Once a value of CL a is obtained, it 
can then be used to determine t} 2 as a function of T]i, (jO x , and from 
equations (B8a) and (BlOa) . Therefore, 


12 


pA \ pA 

+ S SI - S 


(Bll) 


When p is known from experimental measurements , then 


- H ^ < bi2 > 

and oo 2 and r] 2 would be functions only of W], and r) ± , and the values of p 
and CL a - 

The experimental data for a and P are fitted with the tricyclic equa- 
tion (eq. B6) ) by a least-squares procedure using the method of differential 
corrections (ref. 8). To initiate this procedure, starting values for the 
unknown coefficients must be provided. This is accomplished by using a 
modification of Prony f s method (ref. l6) , which will be described next. 

The following theorem is fundamental to Prony*s method: If 

3 

(0 + ia)j = Y, K J eCPjXZ (B13) 

j=i 

where x = xi, x 2 , • • • x n (x locations of equally spaced points. Ax apart) 
and l - 1, 2, • • • n (n = number of points (|3 + ia) ) . Then p + ia satis- 
fies the linear difference equation 


lx 

U 2 = W! - p — 

x y 


consequently. 


i] 2 - 


pA 

31 + 2m CL a 




tOn 


(0 + ia) z+3 + Q 2 (p + ia) i+ 2 + Qi(3 + ia)j +1 + Qo(0 + ia) ^ 
where Q, are constants such that the roots of 

h 3 + Q 2 h 2 + Qih + Qo = 0 


0 

(B13a) 


(B14) 
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are 


hj - 


We have thus minimized 


I 


(0 + ia) z+ 3 + Q2(0 + icc) 7j _ + Qi(0 + ict) Ui + Q 0 (P + ice) 


Z+2 


Z+i 


instead of 


2 


which means simply that the coefficients derived from the procedure are not 
the best possible coefficients but are very good first approximations. To 
obtain values for the coefficients Kj and c pj the following steps would 
normally be taken: Solve the following set of equations by the method of 

least squares for Qi, and Qo- (See eq. (Bl3a).) 


I 


(P + ioo ) 4 + 02(3 + 100)3 + Qi(P + 100)2 + Qo(P + i a )i = 0 

(p + ia)s + Q^( 3 + ico ) 4 + Qi(P + ia )3 + Qo(P + 100)2 = 0 


(0 + ia) + 0 , 2 (. 0 + ice) . + Qi(P + i a )_ _ + Q 0 (P + ^ a )-n „ = 0 

n n-i n -2 u n-3 

With Q 2 y Qi, and Q 0 known , solve the cubic equation 

h 3 + Q 2 h 2 + Qih + Qo = 0 
for roots hi, h 2 , and h 3 where 

hi - 


h s 




I 13 


e 93Ax 
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Therefore, by solving for 9 i, 92 > and cp 3 and then substituting them into 
equation (BI3), one can do a second least-squares operation to obtain Ki, K2, 
and K3. However, since cp 3 is not independent of cp x and 92 for the tricy- 
clic equation in question (i.e., p = f(u)i, w 2 ) ^ g( 9 i, 92) ), a slight 
modification is used. The procedure is: 

Assume initially that 9 3 = 0 and thus h 3 = e° = 1 . Therefore 
( 3 + ia) j = K 1 e' !,lX 2 + K^l * Rj 


and h 3 + Q2I1 2 + Qih + Q 0 = 0 has unity as a root so that 1 + Q2 + Qi + Qo 
If Qo is eliminated between this equation and equation (B 13 a), 


(3 + ia) z+2 - (3 + ia) ? 


Q2 + 


(3 + ia) z+i - (3 + la) l 


Qi 


= ( 3 + ia) l - ( 3 + ia) J+3 

This equation is now solved by least squares for the coefficients Qi and Q2* 
The exponentials and e^ 2 ^ are now found to be the roots of 


h 2 + (Q2 + l)h + (Qi + Q2 +l) — 0 


So , from 


hi 


e 9iAx 


, 92AX 

h 2 = e ^ 


the values of 91 and 92 are found from which 9 3 is calculated; the process 
is repeated with a new 9 3 until the coefficients remain essentially constant 

The coefficient e^j x 2 of Kj in equation (BI3) is tabulated for each 
j since 9j is known. Therefore, Kj can be found by a second least -squares 
procedure . 

The values of a and 3 provided for this routine are either computed 
internally directly from the raw angle data with a table look-up-and- 
interpolation scheme or are input directly with constant increments of 
distance. 
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APPENDIX C 


LIFT EQUATIONS 


If Magnus terms are neglected and if trajectory rather than hody-oriented 
force coefficients are used, Nicolaides 1 equation for lift can be written as 


• • , M 

y + iz 


• • » • • 

y + iz 


- V 2 [c La (P + ia) + ( ol ^ + QLi) (3 + ia) f + (c yo + IQl^ 
+ Cj)(^y “ i/z) 

- V 2 tCL + Cp(7y " i7z)l 


ipt 
e * 


(Cl) 

(Cla) 


where 7y and 7 Z are the projections of the local flight -path angle 7 into 
the y-x and z-x planes, respectively, and z and a data have been corrected 
for gravity effects. Transforming equation (Cl) into x coordinates yields 


y + iz = x(y* + iz') 

. . . . d(y* + iz') 

y + iz = x(y* + iz') + x — 


y + iz = x(y* + iz’) + x 2 (y" + iz") 


Note that 

x s V 

x S KV 2 [ -Cl( 7y ~ i7z) - Cp] , ( K = 
Substituting for x and x into equations (C2) and (C3) yields 

y + iz = v(y* + iz’ ) 

y + iz = KV 2 [-CL(7y - i7z) “ CpUy’ + iz') + V 2 (y" + iz") 


(C2) 


(C3) 


(C2a) 

(C3a) 


If equation (C3a) is substituted into (Cla) with 
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y‘ + iz* s /y - i7 z 


the result is 


y” + iz" = -KC l [1 - (7y - i7 z ) 2 ] 

For small 7, as is the case for ballistic range tests , terms of the order 7 2 
can he neglected. Then 

y" + iz" = -kc l 

y" + iz" = -K £oL a (P + ia) + + Cj,^) (P* + ia*)<i + ^Cy o + iCL 0 ^) e ipX 

(C4) 


After integrating twice with respect to x we obtain 


y + iz = -K 


c y 0 + ±0 L o ) lp p 2 " e ” ) + °L a f Q f Q (P + i^)dx dx 


(d) + Ql^) J (P + iot)dx - (d) ( P Q + ia©)* 


+ (yo + iz o) x + (y Q + iz Q ) 


(C5) 


The integrals are obtained in closed form from the expression for 
P + ice, = Kie'P 1 *- + K 2 ^ 2X + Kse^ 3 ^-. Once the stability routine has been used 
the only unknowns are Cy Q , Cl 0 , (CLq_ + CL&), Jo> z' Q , y Q > and z 0 and 

these appear in a linear fashion and can be found by a straightforward 
least -squares procedure on the experimental y and z data. 
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appendix d 


ROTATION OF COORDINATES 


The purpose of this subroutine is to reduce errors from the small -angle 
assumption inherent in the tricyclic method by rotating the coordinates so 
that the maximum resultant angle of attack falls in either the pitch or yaw 
plane. This is accomplished as follows: After the original gravity-corrected 

angle data have been curve fit (stability subroutine), that solution is exam- 
ined to find the maximum resultant angle of attack that occurs nearest the 
middle of the flight. When it is found, the rotation angle required to bring 
that resultant angle to the yaw plane is determined. 

The experimental 0, ljr, y, and z data are then rotated about the x axis 
through this angle. If the model is rolling rapidly, the peaks of the motion 
will be precessing, and this rotation of coordinates may fail to accomplish 
the stated purpose since successive peaks in the resultant angle of attack 
will not remain in the 3 plane. 

The angle of rotation, 0, and the data in the rotated system are printed 
out for reference. The angle, $, is retained so that when all other parts of 
the data reduction are completed, all information, both data and results, is 
transformed back to the original reference system. 
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APPENDIX E 
SWERVE CORRECTIONS 


The purpose of this subroutine is to correct the flight-path angles a 
and 3 to account for swerve of the model. The relations between the earth - 
fixed angles 9 and t and the flight -path angles a and 3 are 

r? nr 

a = 9 + tan" 1 || (El) 

3 = -t + tan" 1 ^ (E2) 

The procedure for obtaining dz/dx and dy/dx, since they are not known 
a priori, is as follows: First, since they are small, they are assumed equal 

to zero initially and the angle data and swerve data are curve fit (stability 
and lift subroutines). Second, after the initial curve fit to the y and z 
data (lift routine) is obtained, dz/dx and dy/dx are computed at each data 
station and new angles a and 3 are determined. We now recycle through the 
stability and lift routines to obtain new curve fits. From this new curve fit 
to the y and z data, new dz/dx and dy/dx are computed. These are compared 
with the previous set of corrections and if the difference between successive 
corrections is greater than some input factor, say 0.005°, we again recycle 
through the stability and lift subroutines. This continues until all 
corrections are less than 0.005°, at which point the solution is complete . 
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APPENDIX F 


ADDITIONAL FEATURES OF THE PROGRAM 


Dividing Tests Into Segments 

A provision is available for dividing a given set of flight data into as 
many segments of any practical length as desired and treating each segment 
separately. This becomes desirable if the model experiences large velocity 
losses over the length of the range and if the aerodynamics are sensitive to 
changes in Mach number or Reynolds number. It may also be desirable for models 
that have nonlinear aerodynamics with angle of attack and whose maximum ampli- 
tudes change significantly over the length of the range. By analyzing seg- 
ments of such a flight, it is possible to define the variation of aerodynamic 
coefficients with angle of attack from a single test (provided enough data are 
collected to define the motion in each segment). 


Smoothing Procedure 

There is a smoothing procedure that can be used in conjunction with the 
stability and lift routines. Once the iterative scheme involving stability 
and lift has converged (that is, when the angle corrections due to swerve are 
essentially constant) , the differences between the experimental and calculated 
values of a, (3, y, and z are examined at each data station. If any are 
found to be larger than some chosen multiple of the standard deviation of the 
fit in question, then that experimental value is replaced by the calculated 
value and the reduction process for stability and lift is repeated. The 
smoothing process serves two purposes : (l) It tells the experimenter whether 

some data appear to be in error in comparison to the rest; and (2) by repeat- 
ing the calculation it shows him what effect this "bad" data had on the 
derived answers. He can then either check his raw data inputs and resubmit 
the run or conclude that its effect is small enough to ignore. The number of 
times the data are smoothed is determined by the user, usually two or three 
times maximum. If it is done too many times, there is some danger of "walking" 
the least -squares curves in the direction of the replaced values which could 
lead to erroneous results. 


Partially Missing Data 

A provision in the program for weighting the input data (e, i|r, y, z) from 
0 to 1 enables one to place the proper emphasis on each measurement. This is 
particularly useful if, at some data station, part of the measured data are 
either poorly defined or missing (for instance, angular measurement is missing 
in the vertical plane but not in the horizontal) . If the data were missing, a 
weighting factor of 0 would be used (if the data are good, a weighting factor 
of 1 is automatic) . If the data were somewhat uncertain but necessary to 
obtain results, a weighting factor between 0 and 1, say, 0.5, could be used. 
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TABLE I.- TEST DATA 


Test 

M 

RX10" 6 

V, 

ft/ sec 


C L a 


t 

Cm + Cm • 
q_ a 

S +CL a 

deg/ft 

am 

a min 

a rms 







(a) 

Solid models 






[a = 

= 1.25 

in . , m 

£ 0.0229 slug 

x cg 

= 1.95 

% = 

6. 2X10 ~ 4 
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3.59 

-1.259 

-19.1 

-41.5 

12.2 

-.35 

1.80 

.17 

1.28 

1159-1 

2.02 

1.80 

2249 

1.249 

3.36 

-1.490 

-15.8 

-34.0 

18.5 

-.30 

4.79 

.11 

3.55 

1159-2 

1.95 

1.74 

2172 

1.257 

3-35 

-1.426 

-18.6 

-40.9 

33.1 

-.30 

3.86 

.09 

2.74 

1160-1 

1.99 

1.77 

2234 

1.279 

3.31 

-1.495 

-17.5 

-38.3 

14.1 

-.24 

6.84 

.41 

5.00 

1160-2 

1.93 

1.71 

2l6l 

1.284 

3.31 

-1.4i6 

-19.6 

-43.6 

25.8 

-.24 

5.57 

.32 

3.88 







(D) 

Hollow models 
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